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We present a computational method for finding attractors (ergodic sets ol states) of Boolean 
networks under asynchronous update. The approach is based on a systematic removal of state 
transitions to render the state transition graph acyclic. In this reduced state transition graph, all 
attractors are fixed points that can be enumerated with little effort in most instances. This attractor 
set is then extended to the attractor set of the original dynamics. Our numerical tests on standard 
Kauffman networks indicate that the method is efficient in the sense that the total number of state 
vectors visited grows moderately with the number of states contained in attractors. 

PACS numbers: 



I. INTRODUCTION 

Complex disordered systems with many degrees of free- 
dom can often be approximated by Boolean dynamics 
[TJ In particular, gene-regulatory systems in living 
cells [3] have been modeled by Boolean dynamics since 
the 1960s [4]. Each gene is represented by a node with 
two possible states. In this coarse-grained state represen- 
tation, only high (Boolean true = 1) and low (Boolean 
false = 0) chemical concentration of the gene product are 
distinguished. The regulatory biochemical interactions 
between gene product concentrations are captured as log- 
ical rules in the Boolean model. Each unit is assigned a 
Boolean function (truth table) according to which it up- 
dates its state based on the states of other units. In 
recent years, such Boolean models have been shown to 
capture the dynamics of real regulatory systems [SHE] ■ 

The long-term behaviour of Boolean dynamics is of 
particular interest. It is characterized by attractors (er- 
godic sets) as minimal subsets of the state set from which 
the dynamics does not escape. Attractors in a Boolean 
system have been interpreted as distinct cell types in mul- 
ticellular organisms [§1 HO] ■ The computational problem 
of finding all attractors in a Boolean system is difficult. 
Even the simpler problem of deciding whether the sys- 
tem has a fixed point (the smallest possible attractor) is 
NP-complete [TTJ [T^] . In many instances of Boolean net- 
works, however, the state space to be searched may be 
largely reduced [Ml [14] , allowing for attractor detection 
in sparse networks with several dozens of nodes. Known 
methods |15H17| are tailored for dynamics under deter- 
ministic synchronous update. The assumption of fully 
deterministic operation of all nodes in the network may 
not be justified when modeling real regulatory systems. 
In fact, the number and size of attractors change dramat- 
ically when giving up deterministic synchronous update 
[T8] and using stochastic asynchronous update instead 

Here we present a generally applicable exact method 
for attractor search in Boolean dynamics under asyn- 
chronous single-node update. In non-rigorous terms, the 
method works as follows. Departing from the asyn- 



chronous Boolean dynamics to be analyzed, we modify 
the Boolean functions of a subset X of the nodes such 
that they cannot leave a "preferred" state (0 or 1) once 
they have reached it. This reduction of the allowed state 
transitions is guaranteed not to eliminate any of the orig- 
inal attractors: While additional attractors may appear, 
the existing ones may lose some of their states but never 
disappear. In fact, we can choose the set X of nodes 
such that all attractors lose all but one state, i.e. they 
become fixed points. This is the case when each directed 
cycle of the Boolean network contains at least one vertex 
in X. In other words, removal of the nodes in a so- 
called feedback vertex set X leaves the Boolean network 
acyclic. On sparse networks, feedback vertex sets X can 
be found with a cardinality \X\ much smaller than that 
of the whole node set. Then the number of attractors (all 
fixed points) is not greater than 2' X L These attractors 
can be found by systematic enumeration. Now by con- 
struction, each attractor of the original dynamics must 
contain at least one fixed point found for the reduced dy- 
namics. Finally the attractors are found by depth first 
searches seeded at each fixed point of the reduced dy- 
namics. 

The rigorous description of this method in section III 
uses formal notions of Boolean mappings, operators, at- 
tractors and directed graphs. These notions and their 
relations are introduced in section [TTJ Results and con- 
cluding remarks are presented in sections IV and [V] 



II. TECHNICAL BACKGROUND 

A. Boolean dynamics and update operators 

A Boolean dynamical system of n units (or nodes) is 
defined by assigning each node i € {1, . . . ,n} =: n\ a 
Boolean function 



fi : B n -> B 



(1) 



where B = {0, 1} is the set of elementary Boolean states. 
The interaction network underlying this system is ex- 
tracted from the functions fi. There is an edge from 
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node j to node i if and only if function /j explicitly de- 
pends on the j-th coordinate. Put differently, (j, i) is an 
edge if there are state vectors x, y € B" differing only at 
coordinate j such that fi(x) 7^ fi(y). 

The time-discrete dynamics of the system is made pre- 
cise by defining the update mode where synchronous (par- 
allel) update is often used. Then each node computes the 
state at time t + 1 by applying its Boolean function to 
the state vector at time t. A common alternative is the 
fully asynchronous (serial) update mode. At each time 
step, only one node u(t) E {1, . . . , n} is updated while 
all others keep their state. As a generalization, a set of 
nodes potentially containing more than one but not all 
nodes may be chosen individually at each time t. 

We formalize the update modes by defining an update 
operator Ui for each I C {1, . . . , n}. The operator affects 
a Boolean state vector x according to 

(Uix)i — I ^ otherwise ^ 

Then the dynamics is given by the set U of such update 
operators one of which is chosen at each time step for up- 
dating the state vector. For synchronous update, we have 
U = {J7{i,. ,., n }} because the only allowed update involves 
all nodes. Asynchronous update for single nodes is per- 
formed with the operator set U = {U^ : i E {1, . . . , n}}. 
Operators may be applied with different probabilities. 
However, we do not deal with probabilities here because 
the attractors of the system do not depend on them as 
long as each operator in U has positive application prob- 
ability at any time. 



B. Attractors and state transition graph 

A general definition of an attractor is based on the 
state transition graph G — (B n ,T) having all state vec- 
tors B n as its node set. Directed arcs in this graph are 
direct state transitions. So there is an arc (x, y) E T from 
x E B n to y E B n if there is an operator U E U such 
that Ux = y and y 7^ x. The fixed points of the dynamics 
are exactly those state vectors that do not have outgoing 
edges. 

An attractor of the dynamics is a sink component of 
the state transition graph. A non-empty set 5* C B n is a 
sink component of (B n , T) if the following two properties 
are fulfilled. 

(i) For all arcs (x, y) E T with x E S, also y E S. 

(ii) No proper non-empty subset of S has property (i). 

In plain words, S is a minimal non-empty set of state 
vectors from which no arcs point to nodes outside S. By 
A we denote the set of all attractors of the system under 
consideration. Note that A is a set of pairwise disjoint 
sets of state vectors. 



C. Reduced dynamics 

What happens with the set of attractors under small 
modifications of the dynamics? Let us consider the case 
of adding a single arc (x,y) ^ T to the state transition 
graph. Then the modified state transition graph with arc 
set Tl — T U {(x, y)} has an attractor set A' such that 

1. |-4| > \A'\; and 

2. For all Re A 1 there is an 5" E A such that SCR. 

In plain words, the addition of arcs to the state transi- 
tion graph may reduce but not increase the number of 
attractors. Each attractor of the augmented state tran- 
sition graph is contained in an original attractor. Con- 
versely, the removal of arcs may only increase the number 
of attractors; each attractor of the original dynamics is 
contained in an attractor of the dynamics reduced by arc 
removal. This insight suggests to systematically remove 
arcs from the state transition graph to obtain a reduced 
dynamics in which the attractors are easier to find. The 
found attractors can then be used as seeds for the search 
for the attractors in the original dynamics. 

We perform the reduction of the dynamics at the level 
of the update operators. Considering the single-node up- 
dates again, we define the 6-retaining operator for a 
node i and a Boolean state b E {0, 1} by 

U (b) x = { X if Xi = b (3) 
1 y Ux otherwise v ; 

If we replace the operator Ui by the corresponding 6- 
retaining operator, node i can no longer switch from state 
6 to state 1 — 6. This causes the removal of arcs (x, y) 
with Xi = b and t/j = 1—6 from the state transition graph 
while all other arcs remain unaffected. Consequently, ad- 
ditional attractors may be obtained while the existing 
ones become smaller, as discussed in the preceding para- 
graph. 

The goal is now to replace a suitably chosen subset 
of all update operators by state-retaining operators such 
that each attractor shrinks to a single state vector, a fixed 
point. This is certainly the case when the reduced state 
transition graph does not have a directed cycle. Along 
such a cycle, the nodes that change state must do so in 
both directions, flipping from to 1 equally often as from 
1 to 0. On the other hand, these switching nodes induce 
a subnetwork of the Boolean network containing a cycle. 
Then by contraposition, we see how cycles in the reduced 
state transition graph can be suppressed: Each cycle on 
the Boolean network must contain at least one node that 
changes state in only one direction. So we choose a set 
of nodes X that hits each cycle of the graph at least 
once. Such a set is called a feedback vertex set. For each 
node in i E X, the update operator Ui is replaced by a 
state-retaining operator, ensuring that each cycle of the 
Boolean network contains a node that does not flip in 
both directions in the reduced dynamics. 
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FIG. 1: (a) A Boolean network with three nodes and (b) its 
state transition graph G for asynchronous single node update. 
Since G is strongly connected, there is a single attractor com- 
prising all states. The Boolean network may be made acyclic 
by removing node 2 or by removing nodes 1 and 3. Thus 
{2} and {1,3} are the minimal feedback vertex sets (FVS). 
Panel (c) is the reduced state transition graph for FVS {2} 
and 62 = 0, where node 2 cannot leave state X2 = 0. Network 
state 101 is a fixed point of this reduced dynamics. Panel (d) 
shows the reduced state transition graph for the larger FVS 
{1,3} with retained states 61 = 63 = 0. Here the reduced 
dynamics has a set of three fixed points, P — {000, 001, 100}, 
inside the original attractor. 

Since X is a feedback vertex set of the Boolean net- 
work, the subnetwork induced by the remaining nodes 
Y = rtj \ X does not have cycles. It is a feed-forward 
network. Therefore there is a topological sorting of Y: 
The indices of the nodes in Y can be permuted such that 
arcs between nodes in Y go only from a lower index to a 
higher index. Given that the system is in a fixed point, 
the states of the nodes in X are sufficient to calculate the 
states of the remaining nodes Y as well. This is useful 
for finding all fixed points of the reduced dynamics: For 
each state vector on X, complete it to a state vector on 
the whole network and check if this state vector is a fixed 
point or not. 

The set P of all fixed points of the reduced dynamics 
serves as a starting point for constructing the attractor 
set A of the original dynamics. Each attractor S G A 
contains at least one of the state vectors in P. Thus 
for each x* G P we calculate the set of all state vectors 
reachable from the x*. If this set contains another ele- 
ment of P then x* may be discarded. Otherwise this set 
is an attractor. 



III. METHOD FOR FINDING ATTRACTORS 

The computational method for finding attractors falls 
into three stages: (A) Establish a feedback vertex set that 
defines which update operators are made state-retaining. 



(B) Find the fixed point set A* of the reduced dynamics 
by enumeration (C) Traverse the original state transition 
graph departing from the state vectors in A*. A small 
example is illustrated in Figure [T] 

A. Feedback vertex set 

A feedback vertex set X of the given Boolean network 
is determined as follows. Initialize X as the whole node 
set n\ . Loop: Draw a node i G X at random (flat distri- 
bution). If X \ {{] is a feedback vertex set, set X \{i}. 
Repeat this loop until X does not have a proper subset 
being a feedback vertex set. This very simple method 
tends to generate a small but not necessarily globally 
minimal feedback vertex set. 



B. Fixed points of reduced dynamics 

In addition to the feedback vertex set X, the retained 
state hi needs to be determined for each i G X. Here we 
choose bi to be the state that /j assumes for most values 
of the argument. If there is a tie between 0s and Is, bi is 
drawn at random from B with equal probabilities. 

Without loss of generality we assume an indexing of the 
nodes such that X = {1, . . . , m} and m + 1, m+ 2, . . . , n 
is a topological sorting of the feed-forward subnetwork 
induced by the remaining nodes n\ \ X =: Y. After 
initializing the set P of fixed points as the empty set, we 
perform the following nested loops. The outer loop runs 
over all partial state vectors x G B m . The first inner loop 
runs over the vertices in i G Y in the order of topological 
sorting, calculating the state Xi = fi{x\, . . . , Xi-i). After 
finishing the first inner loop, a second inner loop checks 
if 

C/f i] x = x (4) 

for all i G X. If yes, a; is a fixed point of the reduced 
dynamics so x is included in P. Otherwise x is discarded. 



C. Original attractors 

After initializing the set of attractors A as the empty 
set, attractor finding is performed as the following steps. 
(1) Draw an element x* G P and remove it from P. (2) 
Perform a depth first search on the state transition graph 
starting at x* G P. (3) If the search encounters an el- 
ement x G P, it terminates immediately and the search 
result is discarded. Otherwise the search runs until no 
further unvisited state vectors are found. Then the set S 
of state vectors visited during the search is included in A 
as an attractor. (4) If P is not empty, resume at (1), oth- 
erwise all attractors have been found. Figure [T] shows a 
case where the reduced dynamics has several fixed points 
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FIG. 2: Computational effort of attractor finding measured 
as the ratio r between the number of states t calculated and 
the number of states s actually contained in attractors. The 
inset shows a linear-log plot of the same data pointing out the 
weak increase of r with s. The method is applied to Random 
Boolean networks with K — 2 inputs per node (critical Kauff- 
man networks). Ensembles contain 10 4 independently gener- 
ated network instances for each system size n G {13,25,63}. 
A realization yields a data point (s, r) with s and r as defined 
in the main text. Each plotted point is an average over data 
points with s G [2 k , 2 k+1 - 1], k = 0, 1, 2, . . . (logarithmic bin- 
ning). For n = 63, only 8551 data points are used. On the 
remaining instances, the computation runs out of memory. 
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FIG. 3: Scatter plot of computational steps t versus number 
of states on attractors s. Each data point represents a network 
realization with n = 63 nodes. See caption of [2] for details. 



IV. RESULTS 

Standard random Boolean networks with K — 2 inputs 
per node are generated as test instances for the method 
as follows. Each node i is assigned randomly and in- 
dependently one out of the 16 Boolean functions of two 
variables and a pair of nodes from which node i receives 
input. 

We test the method with various system sizes up to n — 
63 and several independent realizations of a network. As 
a first qualitative result we find that the computational 
time of steps A and B of the method (finding X and P) 
is negligible compared with the final effort of searching 
the state transition graph. Therefore we measure the 
computational cost in terms of the number t of successor 
states computed. Note that t may be larger than 2™ 
because a state with several predecessor states may be 
computed more than once. A lower bound for t is the 
number s of state vectors actually contained in attractors 
where 



seA 



(5) 



contained in the same original attractor so search results 
would be discarded in step (3). 



The method is efficient if the ratio r := t/s is small. 
Figure [2] shows that low values of r are obtained for in- 
termediate values of s. With increasing s, the value of 
r grows moderately. The functional form might be loga- 
rithmic, r ~ logs. More extensive numerical simulations 
or analytical estimates are required to clarify the growth. 
The scatter plot in Figure [3] shows that fluctuations of t 
become small for large s. 

At network sizes n > 30, the computation does not 
succeed for all the 10000 instances because it exceeds 
the allocated memory. At n = 63, a fraction of 14.5% 
is unsuccessful because the depth first search of G does 
not stay within a maximum depth of 3 x 10 7 . Restarting 
with a different random number for the 1449 unsuccessful 
instances leads to a complete computation in merely 19 
additional cases. This suggests that most of the failing 
instances are intrinsically difficult and do not fail because 
of an unfortunate choice of the feedback vertex set. 



V. CONCLUDING REMARKS 

As a proof of concept we have introduced and used 
the method in its simplest form. Several extensions of 
the method may help to increase the efficiency and allow 
computation of attractors for larger and denser networks. 

Smaller feedback vertex sets may be found by replac- 
ing our simple greedy approach with a more advanced 
method. Having smaller sets X tends to reduce the set 
of fixed points P of the reduced dynamics. It remains 
to be seen if smaller P leads to systematically shorter 
searches of the state transition graph in the last phase of 
the computation. 

One could attempt to explicitly guide depth first search 
towards other elements in P such that futile searches 
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would terminate faster. The search then should first 
choose as the next state vector a predecessor that reduces 
the distance to another element in P. Alternatively, sys- 
tematically different graph traversals, e.g. breadth first 
search, might be tried out. 

As the greatest limitation of the method, we experi- 
enced the need to store all the state vectors of an at- 
tractor during the traversal of the state transition graph. 
Therefore for some of the realizations at large n the com- 
putation ran out of memory. Space efficiency might be 
gained by a compression of sets of state vectors during 
the traversal. Here it may help that attractors of sparse 
networks often have a combinatorial product structure 

Can the method be modified to work with syn- 



chronously updated Boolean systems as well? Any syn- 
chronous system can be emulated with asynchronous up- 
date under sufficient extension of the state space [21]. It 
would be interesting if the pruning of arcs in the state 
transition graph is a general principle applicable to a 
larger class of discrete dynamics. 
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